********************************************************************************
*Estimating survival functions for Beardsley's original model
********************************************************************************


use book-basedata-replication, clear

*splitting at failure time
stsplit, at(failures)

*recalculate the interaction with time
foreach var in med2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 {
replace `var'_t= `var'*_t
}

*Beardsley's original model
stcox med2 med2_t prevcris2 viol2 crisdur2 jointdem victory2 contig2 /* 
*/		if _t<3650,  strata(order5) cluster(dyadno) efron nohr
matrix b=e(b)

*predict baseline hazard contribution
predict baseline, basehc

*calculate mean of each variable
foreach var of varlist med2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 {
quietly sum `var' if e(sample)
local s_`var'=r(mean)
}

*aggregate on analysis time level and create baseline hazard for each stratum
preserve
collapse (mean) baseline, by(_t order5)
reshape wide baseline, i(_t) j(order5)

*for each stratum...
forvalues t=1/5 {
*...calculate survival function if no mediation and mean covariate values...
gen base_h_mean`t'=baseline`t'*exp(b[1,3]*`s_prevcris2' +b[1,4]*`s_viol2' + /* 
*/		b[1,5]*`s_crisdur2' +b[1,6]*`s_jointdem' +b[1,7]*`s_victory2' + /* 
*/		b[1,8]*`s_contig2') if _t<547
replace base_h_mean`t'=baseline`t'*exp(b[1,3]*`s_prevcris2' +b[1,4]*`s_viol2' + /* 
*/		b[1,5]*`s_crisdur2' +b[1,6]*`s_jointdem' +b[1,7]*0 + /* 
*/		b[1,8]*`s_contig2') if _t>=547
gen base_H_mean`t'=sum(base_h_mean`t')
gen base_S_mean_B`t'=exp(-base_H_mean`t')

gen med_h_mean`t'=baseline`t'*exp( /*b[1,1]+b[1,2]*_t+ <- no mediated outcome before day 547
*/		b[1,3]*`s_prevcris2' +b[1,4]*`s_viol2' + /* 
*/		b[1,5]*`s_crisdur2' +b[1,6]*`s_jointdem' +b[1,7]*`s_victory2' + /* 
*/		b[1,8]*`s_contig2') if _t<547
replace med_h_mean`t'=baseline`t'*exp(b[1,1]+b[1,2]*_t+ /*
*/		b[1,3]*`s_prevcris2' +b[1,4]*`s_viol2' + /* 
*/		b[1,5]*`s_crisdur2' +b[1,6]*`s_jointdem' +b[1,7]*0 + /* 
*/		b[1,8]*`s_contig2') if _t>=547
gen med_H_mean`t'=sum(med_h_mean`t')
gen med_S_mean_B`t'=exp(-med_H_mean`t')
}

*keep survival estimates and time variable and store
keep base_S_mean_B* med_S_mean_B* _t
rename _t _t_new2
save my_surv_curve1, replace

*merge with analysis data
restore
capture drop _merge
merge using my_surv_curve1

label var base_S_mean_B1 "No mediation"
label var med_S_mean_B1 "Hypothetical mediation 1.5 years after crisis ended"

*graph results for each stratum
twoway line base_S_mean_B1 med_S_mean_B1 _t_new2 if _t_new<3650, sort c(J J) /* 
*/		lcol(gs7 gs10) saving(order1_dem_change, replace) legend(col(1) region(lwidth(none)))/* 
*/		t2title("1 previous crisis") xtitle("Days at peace") ylab(.4(.2)1) ytitle("")
forvalues t=2/4 {
twoway line base_S_mean_B`t' med_S_mean_B`t' _t_new2 if _t_new<3650, sort c(J J) /* 
*/		lcol(gs7 gs10) saving(order`t'_dem_change, replace) legend(off)  /* 
*/		t2title("`t' previous crises") xtitle("Days at peace") ylab(.4(.2)1) ytitle("")
}
twoway line base_S_mean_B5 med_S_mean_B5 _t_new2 if _t_new<3650, sort c(J J) /* 
*/		lcol(gs7 gs10) saving(order5_dem_change, replace) legend(off)  /* 
*/		t2title("5 or more previous crises") xtitle("Days at peace") ylab(.4(.2)1) ytitle("")


grc1leg order1_dem_change.gph order2_dem_change.gph order3_dem_change.gph /* 
*/		order4_dem_change.gph order5_dem_change.gph, ycommon /* 
*/		l2("Probability of survival") /* 
*/		leg(order1_dem_change.gph) 
graph export Beardsley_delayed_survival_atmean_strata.pdf, as(pdf) replace

